*This code must be run within the main do file. It will not run correctly on its own. 



	use "$data/Replication_data.dta", replace
	
	local a_dry = " prec_dry below40_dry FND_Max40to45_dry FND_Max45to50_dry FND_Max50to55_dry FND_Max55to60_dry FND_Max60to65_dry FND_Max75to80_dry  FND_Max65to70_dry "
	
	local a_wet = " prec below40_wet FND_Max40to45_wet FND_Max45to50_wet FND_Max50to55_wet FND_Max55to60_wet FND_Max60to65_wet FND_Max75to80_wet  FND_Max65to70_wet "
	
	local b_dry = " c.Dividends#c.prec prec    c.Dividends#c.below40_dry c.Dividends#c.FND_Max40to45_dry c.Dividends#c.FND_Max45to50_dry c.Dividends#c.FND_Max50to55_dry c.Dividends#c.FND_Max55to60_dry c.Dividends#c.FND_Max60to65_dry c.Dividends#c.FND_Max75to80_dry  below40_dry FND_Max40to45_dry FND_Max45to50_dry FND_Max50to55_dry FND_Max55to60_dry FND_Max60to65_dry FND_Max75to80_dry FND_Max65to70_dry c.Dividends#c.FND_Max65to70_dry"	
	
	
	local b_wet = " c.Dividends#c.prec prec    c.Dividends#c.below40_wet c.Dividends#c.FND_Max40to45_wet c.Dividends#c.FND_Max45to50_wet c.Dividends#c.FND_Max50to55_wet c.Dividends#c.FND_Max55to60_wet c.Dividends#c.FND_Max60to65_wet c.Dividends#c.FND_Max75to80_wet  below40_wet FND_Max40to45_wet FND_Max45to50_wet FND_Max50to55_wet FND_Max55to60_wet FND_Max60to65_wet FND_Max75to80_wet  FND_Max65to70_wet c.Dividends#c.FND_Max65to70_wet"	
	
	foreach n in mean_Div4 { 
		cap rename Dividends x
		rename `n' Dividends
		label var lny_nonag "`n'"
		
			reghdfe lny_nonag days90plus_dry FND_Max85to90_dry FND_Max80to85_dry  c.Dividends#c.days90plus_dry c.Dividends#c.FND_Max85to90_dry c.Dividends#c.FND_Max80to85_dry Dividends`b_dry', absorb(i.fips i.year i.state#c.year) cluster(i.fips#i.year) 
				
		
		bysort DivCentile: egen Divmean=mean(Dividends)
			levelsof Divmean, local(valuesa)
			quietly margins, dydx(days90plus_dry) at(Dividends=(`valuesa'))
				mat a1=r(table)
				mat b1=r(b)
					g Out=.
					g UL=.
					g LL=.
		
			forvalues q=1(1)100{
				replace Out=b1[1,`q'] if DivCentile==`q'
			}
		
			forvalues q=1(1)100{
				replace UL=a1[6,`q'] if DivCentile==`q'
			}

			forvalues q=1(1)100{
				replace LL=a1[5,`q'] if DivCentile==`q'
			}
			
			twoway  (rcap LL UL DivCentile, msize(vtiny) lc(gs13) scheme(s1color) ///
				yline(0) ytitle("Change in annnual payroll per 90F day (%)") title("") ///
				xtitle("Centile of per capita dividend payments") ylabel(-0.001 "-0.1" -0.0005 "-0.05" 0 "0" 0.0005 "0.05" 0.001 "0.1" 0.0015 "0.15") xsc(reverse) legend(off)) (scatter Out DivCentile)
				
			
							
							graph export "$output/images/Fig1a.png", as(png) replace 				
				
		
		drop Out UL LL Divmean
		
		bysort PovCentile: egen PovMean=mean(Dividends)
			levelsof PovMean, local(valuesb)
			quietly margins, dydx(days90plus_dry) at(Dividends=(`valuesb'))
				mat a2=r(table)
				mat b2=r(b)
					g Out=.
					g UL=.
					g LL=.
			forvalues q=1(1)100{
				local qq=(`q'-101)*-1
				replace Out=b2[1,`qq'] if PovCentile==`q'
			}
		
			forvalues q=1(1)100{
				local qq=(`q'-101)*-1
				replace UL=a2[6,`qq'] if PovCentile==`q'
			}

			forvalues q=1(1)100{
				local qq=(`q'-101)*-1
				replace LL=a2[5,`qq'] if PovCentile==`q'
			}
			
			twoway  (rcap LL UL PovCentile, msize(vtiny) lc(gs13) scheme(s1color) ///
				yline(0) ytitle("Change in annnual payroll per 90F day (%)") title("") ///
				xtitle("Centile of population below 200% of 2010 poverty line") ylabel(-0.001 "-0.1" -0.0005 "-0.05" 0 "0" 0.0005 "0.05") legend(off)) (scatter Out PovCentile)
						
							graph export "$output/images/Fig1b.png", as(png) replace 				
						
			
			drop Out UL LL
			
			rename Dividends `n'
	}
	
